Nonlinear waves in a cylindrical Bose-Einstein condensate 
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We present a complete calculation of solitary waves propagating in a steady state with constant 
velocity v along a cigar-shaped Bose-Einstein trap approximated as infinitely-long cylindrical. For 
sufficiently weak couplings (densities) the main features of the calculated solitons could be captured 
by effective one-dimensional (ID) models. However, for stronger couplings of practical interest, 
the relevant solitary waves are found to be hybrids of quasi- fD solitons and 3D vortex rings. An 
interesting hierarchy of vortex rings occurs as the effective coupling constant is increased through 
a sequence of critical values. The energy-momentum dispersion of the above structures is shown to 
exhibit characteristics similar to a mode proposed sometime ago by Lieb within a strictly ID model, 
as well as some rotonlike features. 

PACS numbers: 05.30.Jp, 03.75.Fi, 05.45.Yv 



I. INTRODUCTION 

Solitary waves that may occur in a Bose-Einstein Con- 
densate (BEC) have been traditionally discussed in terms 
of the classical Gross-Pitaevskii (GP) model which is ap- 
propriate for the description of weakly correlated sys- 
tems 0. For instance, a simple soliton was obtained by 
Tsuzuki H in a homogeneous ID model, while Zakharov 
and Shabat f| developed inverse-scattering techniques 
for the study of multisolitons. Interestingly, the elemen- 
tary soliton proved to be relevant for an accurate semi- 
classical description El |Bj of an intriguing mode proposed 
earlier by Lieb Q in a full quantum treatment of a ID 
Bose gas based on the Bethe Ansatz 0. 

The above developments had long remained purely the- 
oretical because of the absence of a physical realization 
of a strictly ID Bose gas. Nevertheless, the picture has 
significantly changed with the recent observation of sim- 
ilar coherent structures in confined BECs of alkali-metal 
atoms ^ |j| . The very method of experimental produc- 
tion of solitary waves (phase imprinting) was inspired by 
the analytical structure of the ID soliton, while various 
effective ID models have been developed for their theo- 
retical investigation @ [□], g f0| [b§ [0§. On the 
other hand, the actual stability of the theoretically pre- 
dicted ID solitary waves should be questioned within the 
proper 3D environment of realistic traps |12, p], An 
important step in that direction was the experimental 
observation jl7| that a dark soliton initially created in a 
finite trap eventually decays into vortex rings, as is also 
predicted by a numerical solution of the corresponding 
initial- value problem in a 3D classical GP model [jl3) . 

Therefore, it is important to carry out a calculation of 
potential nonlinear modes without a priori assumptions 
about their effective dimensionality. One could envisage 
a picture in which the actual solitary waves are hybrids 
of quasi- ID solitons and 3D vortex rings. It is the aim 
of the present paper to make the above claim precise by 
calculating solitary waves that propagate along a cylin- 
drical trap in a steady state with constant velocity v. 



Our approach was motivated by the calculation of vortex 
rings in a homogeneous BEC due to Jones and Roberts 
fL8[ and a similar calculation of semitopological solitons 
in planar ferromagnets fl9|| . 

We have already described the main result of this work 
in a recent short communication ]2Q| but a substantial 
elaboration is necessary in order to appreciate its full 
significance. Thus the problem is formulated in Sec. II 
where we also present a brief but complete recalculation 
of the ground state and the corresponding linear (Bogoli- 
ubov) modes for comparison. A detailed calculation of 
nonlinear modes is given in Sec. Ill and the main conclu- 
sions are summarized in Sec. IV. 



II. FORMULATION AND LINEAR MODES 

The physical picture that we have in mind is a slight 
idealization of the experiment in Ref. Thus we 

consider a cigar-shaped trap filled with atoms of mass 
to. The transverse confinement frequency is denoted 
by uj± and the corresponding oscillator length by a±_ = 
(h/mui^) 1 / 2 . The longitudinal confinement frequency 
u>\i is assumed to be much smaller than uj±, hence we 
make the approximation of an infinitely-long cylindrical 
trap with u>» —0. Accordingly, complete specification of 
the system requires as input the average linear density 
v which is the number of atoms per unit length of the 
cylindrical trap. Finally, we consider the two dimension- 
less combinations of parameters: 
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7-L = va ±- 



(1) 



where a is the scattering length related to the coupling 
constant as usual by Uo = ATrh 2 a/m. 

Now, in the actual experiment of Ref. || , the trap is 
filled with 87 Rb atoms, the transverse frequency is chosen 
as lj± = 2ir x 425 Hz, the oscillator length is calculated 
to be a± « 0.5/im, and the estimated linear density lies 
in the range v> 0.1 atoms/A. Also taking into account a 
scattering length a sa 50 A, the dimensionless parameters 
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(|l|) take values in 7 > 5 and j± > 5 x 10 2 . In fact, our 
subsequent calculations will be carried out for a much 
wider range of the above parameters. Therefore, apart 
from the idealization of a cylindrical trap, our results are 
fairly realistic and could be applied to a number of cases 
of experimental interest. 

It is useful to introduce rationalized units through the 
rescalings 



t 



,1/2 



(2) 



The energy functional extended to include a chemical 
potential is then given by 

W = - J [V** V* + p 2 + g(***) 2 - 2/j, dV, 

(3) 

where g = Attj and p 2 = x 2 +y 2 . Eq. (m yields energy W 
in units of r y±{hu;±) whereas the chemical potential /x is 
measured in units of Hlo±. The corresponding rational- 
ized equation of motion reads 



1 1 , 



(4) 



and depends only on the dimensionlcss coupling constant 
7, because g = Attj and the chemical potential p = fJt('j) 
is fixed by the requirement that the system carry in its 
ground state a definite average linear density v. 

An important first step is thus to obtain accurate infor- 
mation about the ground-state wave function , E r = 1 E r o(p) 
which is normalized according to 



2tt P dp |^o| 2 = 1 



(5) 



to conform with our choice of rationalized units. The 
wave function ^o(p) is numerically calculated as the min- 
imum of the energy functional, under the constraint (|^) 
that fixes the chemical potential ^(7), by a variant of a 
relaxation algorithm plj . Explicit results are illustrated 
111 Fig. I for some typical values of 7 where we also quote 
the corresponding values of the chemical potential. 

The preceding numerical determination of the ground 
state will provide the basis for all subsequent calcula- 
tions. However, it is worth mentioning here some limit- 
ing cases where the ground state is known analytically. 
At 7 = 0, 



= -p72 



(6) 



and the chemical potential degenerates to p = 1 . In the 
opposite limit, 7 S> 1, one may use the Thomas-Fermi 
(TF) approximation p2| 
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(7) 



for < p < R ± = 27 1 / 4 , and * = for p > R ± . The 
chemical potential is given accordingly by p^2^/ 1 / 2 . A 




FIG. 1: Radial dependence of the ground- state wave function 
for four representative values of the dimensionless coupling 
constant 7 = 0, 1, 10 and 100. The corresponding values of the 
chemical potential were calculated to be p = 1, 2.2571, 6.4324, 
and 20.0431, in units of hu±. The inset compares the accurate 
numerical solution (solid line) with the TF approximation (Q) 
(dotted line) for the strong coupling 7 = 100. Distance is 
measured in units of a±. 



comparison with the accurate numerical solution is shown 
in the inset of Fig. [j] for 7 = 100. In fact, as we shall 
see shortly, the TF approximation provides a reasonable 
description of some quantities of physical interest even 
for 7 ~ 1 . 

We will also need some information from the lin- 
ear (Bogoliubov) modes which have already been calcu- 
lated in the literature to varying degree of completeness 
p3| , p4| , |25| . Here we employ a numerical algorithm of 
our own briefly described as follows. Equation is lin- 
earized by inserting ^ = x &o + f + i\ where ^o(p) is the 
calculated ground-state wave function while / — f(r,i) 
and x = x( r > t) are real functions that account for small 
fluctuations around the ground state. It is somewhat 
more convenient to use the linear combinations a = f + x 
and b = f — x which satisfy the linearized equations 



d_ 

dt 



= M 



where 



1 



M = 



1 



-V -D 
D V 



V = g^i D = --A + -p 2 + 2V-p 



(8) 



(9) 



Our task is then to calculate the spectrum of the differen- 
tial operator M whose eigenvalues are purely imaginary 
and come in pairs ±iu where uj is the sought after phys- 
ical frequency. 
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FIG. 2: The lowest branch in the Bogoliubov spectrum for 
four representative values of the dimensionless coupling con- 
stant 7 = 0, 1, 10 and 100. The frequency U) is measured in 
units of ui± and the wave number q in units of l/a±. The 
corresponding values of the speed of sound were calculated to 
be c = 0, 0.95, 1.77, and 3.17, in units of a±u±. 



FIG. 3: Speed of sound c in units of a±u)± as a function of 
the dimensionless coupling constant 7. Open circles stand 
for our numerical data, the solid line for the TF asymptote of 
Eq. ( fll(|. a nd the dotted line for the weak-coupling asymptote 

ofEtnyj. 



We restrict attention to axially-symmetric waves that 
propagate along the z axis with wave number q. The 
Laplace operator is then replaced by 



A 



dp 2 



ld_ 

pdp 



(10) 



and the amplitudes a and b may be assumed to depend 
only on the radial distance p. A finite-matrix approx- 
imation of the operator M is obtained by expanding 
both a and b in terms of a basis set of non-orthogonal 
Gaussian wave packets with randomly chosen oscillator 
lengths |^6| . It is also prudent to enlarge the basis set 
by including the ground-state wave function x I'o(/c) it- 
self, in order to directly account for the zero (Goldstone) 
mode associated with the number symmetry. The result- 
ing algorithm is then quite efficient and provides stable 
approximations of the low-lying eigenvalues even if we 
include a small number of basis elements. 

In Fig. ||we present explicit results for the lowest eigen- 
frequency uj = tu(q) for the same set of coupling constants 
as in Fig. |l|. At 7 = 0, u(q) reduces to the free-particle 
quadratic dispersion u — q 2 /2, as expected. At nonzero 
7, the dispersion becomes linear near the origin, cj»c|g|, 
where c is the speed of sound for which explicit values 
are also quoted in Fig. ||. Finally, we note that our re- 
sults are in apparent agreement with the Bogoliubov dis- 
persion calculated earlier within the TF approximation 
pM 24 as well as numerically p3 - even though a dif- 



the latter reference. 

The speed of sound is a quantity of special physical 
interest and will also play an important role in the theo- 
retical development of Sec. III. Hence we have carried out 
a calculation for a wider set of coupling constants and the 
results are summarized in Fig. [j[ It is interesting that 
our accurate numerical results arc consistent with the TF 
approximation §1 |J, |tJ 



vV4 



(11) 



even for values of 7 as low as 1, where the error is about 
5%, whereas the error is reduced to less than 1% for 
7^10. This fact is especially important because Eq. (O) 
was employed for the analysis of experimental data p8j . 
The relative accuracy of this approximation progressively 
deteriorates in the region 7 < 1, but a new asymptote, 
namely 



(27) 



1/2 



(12) 



ferent parameterization of the spectrum was employed in 



was predicted to be reached for sufficiently weak cou- 
plings Q . The weak-coupling approximation (|l2|) is ac- 
tually consistent with our numerical data for 7 < 1/4, as 
is shown in the inset of Fig. ||. However, we should add 
that the linear part of the Bogoliubov dispersion becomes 
very narrow in this region of couplings. 
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III. NONLINEAR WAVES 

We now turn to the calculation of axially-symmetric 
solitary waves traveling along the z axis in a steady state 
with constant velocity v. These are described by a wave 
function of the form <3> = ^(p, £), with £ = z — vt, which 
is inserted in Eq. (Q) to yield the stationary differential 
equation 



9* 
A 



2 2 H 



■ff(**¥)*-/i¥,(13) 



<9p 2 



1_9 

P dp 



The wave function must vanish in the limit p — > oo, 
thanks to the transverse confinement, while the condi- 
tion 



km 1*09,01 = 1*0001 



(14) 



enforces the requirement that the local particle density 
coincide asymptotically with that of the ground state cal- 
culated in Sec. II. But the phase of the wave function is 
not fixed a priori at spatial infinity except for a mild 
restriction implied by the von Neumann boundary con- 
dition 



9^ „ 
hm — = 



(15) 



adopted in our numerical calculation. Our task is then to 
find concrete solutions of Eq. (|l^) that satisfy the bound- 
ary conditions just described. 

An important check of the numerical calculation is pro- 
vided by the virial relation 

f Tl <9** d^> a 
vV = J -^^+p 2 *** + |(***) 2 ~M*** dV, 

_ (16) 

obtained by standard scaling arguments fl9|| . Here V is 
the linear momentum given by the usual definition 



V 



1 

2i 



dz 



dz 



* \ dV = 



n^dV, (17) 
dz 



and is measured in units of hv = 7j_(7i/ci_l). In the sec- 
ond step of Eq. ( [j"7|) we employ hydrodynamic variables 
defined from 



<n e 



(18) 



where n= |*| 2 is the local particle density and the phase 
(f> may be used to construct the velocity field u = V0. 

Numerical solutions of Eq. (|l^) are obtained by an 
iterative Newton-Raphson algorithm |ll| |l9| briefly de- 
scribed as follows. Suppose that * = ^- lrL is an initial 
rough guess for the solution at some velocity v. We 
then insert in Eq. ( |l3| ) the configuration * out = \&j n + X 



and keep terms that are at most linear in the ampli- 
tude X. Thus we derive an inhomogeneous differential 
equation of the form LX = Y where the linear opera- 
tor L and the source Y are both calculated in terms 
of $! m . We solve this linear system for X = L~ l Y to 
obtain * ut = + L~ X Y which is used as input for 
the next iteration until convergence is achieved at some 
specified level of accuracy. The procedure is repeated 
by incrementing the velocity to a different value, typi- 
cally in steps of <5u = ±0.01, using as input the converged 
configuration obtained at the preceding value of the ve- 
locity. Therefore, the main numerical burden consists 
of constructing a finite-matrix lattice approximation of 
the linear operator L which is then inverted by standard 
routines appropriate for sparse linear systems. 

The Newton-Raphson algorithm typically converges af- 
ter a few iterations and the final configuration is indepen- 
dent of the specific choice "fin- But it is also clear that 
the algorithm will not converge for most choices of ^i n . 
Hence it is important to invoke an educated guess for the 
input configuration provided by the product Ansatz: 



* in = [ci - i c 2 tanh(c 3 £)] * (p) 



(19) 



which capitalizes on the analytically known solitary wave 
in the homogeneous ID model || and the ground-state 
configuration v I / o(p) numerically calculated in Sec. II. The 
constants c\ , c 2 and C3 are definite functions of the ve- 
locity v within the strictly ID model, but such precise 
relations need not be invoked for our current purposes 
except for the normalization condition c 2 + c 2 — 1 that 



is necessary to enforce the boundary condition (|14|). In 
other words, the above constants are treated here as trial 
parameters until we achieve convergence for a specific ve- 
locity v. A corrolary of the preceding discussion is that 
the converged configuration does not depend on the pre- 
cise choice of those parameters, and it is certainly not 
in the form of a product Ansatz often employed for the 
derivation of effective ID models pi Finally, we 

note that the Ansatz (|l^) satisfies the parity relations 

Re*(p.O = Re*(>,-£), lm*(p,0 = -Im*(p,-^) 

_ (20) 
which are compatible with Eq. (O) and are actually sat- 
isfied by all solutions constructed in the present paper. 

We begin with the special case of the relatively weak 
coupling 7=1 for which the speed of sound was calcu- 
lated to be c = 0.95 in Sec. II. The simplest possibility is 
to first attempt to derive a static (v = 0) soliton starting 
with the input configuration (|l9| ) applied for, say, c\ = 
and C2 = 1 = C3. Indeed, the algorithm quickly converges 
to a wave function with a nontrivial imaginary part but 
vanishing real part. The velocity is then incremented to 
positive values in steps of 6v = 0.01 and the correspond- 
ing wave functions acquire also a nontrivial real part. 
The process may be continued until the velocity v ap- 
proaches the speed of sound c beyond which the solitary 
wave ceases to exist. An equivalent sequence of solitary 
waves with velocities in the range — c<v< is obtained 
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FIG. 4: Solitary wave for 7 = 1 and two representative values 
of the velocity: v~c/2 (I) andu = (II). We display the radial 
dependence of the real and the imaginary part of the wave 
function for various positive values of £ in steps of 5^ = 0.1. 
The corresponding results for negative £ are obtained through 
the parity relations (M) . 




FIG. 5: Contour levels of the local particle density n= | \I/ 1 2 for 
7 = 1, in a plane that contains the z axis and cuts across the 
cylindrical trap. The complete 3D picture may be envisaged 
by simple revolution around the z axis. Regions with high 
particle density are bright while regions with zero density are 
black. The two special cases considered are the same as in 
Fig. 0. 



either by starting again with the v = soliton and pushing 
it to negative velocities or, simply, by taking the complex 
conjugate of the wave function calculated for < v < c, 
since 



-v, * — * 



(21) 



is an obvious symmetry of Eq. (|i~3|). A detailed illus- 
tration of the calculated solitary wave function is given 
in Fig. [| for two representative values of the velocity: 
v~c/2 and v — 0. 

A partial but more transparent illustration is given in 
Fig. H which depicts the level contours of the local parti- 
cle density n= l^l 2 for the two special cases considered 
in Fig. |]. In words, the calculated solitary wave is a 
mild soundlike disturbance of the ground state when \v\ 
approaches the speed of sound c, while it becomes an in- 
creasingly dark soliton with decreasing \ v\ and reduces to 
a completely dark (black) soliton at v = 0. 

It is now important to calculate the energy-momentum 



dispersion of the solitary wave. The excitation energy is 
defined as 



E = W-W ( 







(22) 



where both W and Wq are calculated from Eq. ap- 
plied for the solitary wave ^(p, £) and the ground state 
^o(p), respectively. The presence of the chemical poten- 
tial in Eq. (||) provides the compensation that is neces- 
sary to compare energies of states with the same number 
of particles. Similarly, the relevant physical momentum 
is not the linear momentum V of Eq. ( |l7| ) but the im- 
pulse Q defined in a manner analogous to the case of a 
homogeneous gas m, H , 



Q = I (n-n )-^dV = V 



(23) 



2irpdp n {p)[(j){p, z = oo) — (f)(p, z = — 00)] 



where no = | \I/o (p) 1 2 is the ground-state particle density 
and 5<f> is now the weighted average of the phase differ- 
ence between the two ends of the trap. The delicate dis- 
tinction between linear momentum and impulse has been 
the subject of discussion in practically all treatments of 
classical fluid dynamics ^| and continues to play 
an important role in the dynamics of superfluids [Q. 
Here we simply postulate the validity of the definition 
of impulse in Eq. ( p3| ) and note that the corresponding 
group- velocity relation 



dE 
~dQ 



(24) 



is satisfied to an excellent accuracy in our numerical cal- 
culation and thus provides a highly nontrivial check of 
consistency. In turn, the virial relation (|l^) is verified 
using the standard definition of the linear momentum V 
in Eq. (|lj), as expected. We finally note that the same 
phase difference Scj) which is important for experimental 
production of solitary waves through phase imprinting 
H ^| is also crucial for the calculation of the impulse. 

The dispersion E — E(Q) calculated for the complete 
sequence of solitary waves with velocities in the range 
—c<v<c is illustrated in Fig. ^. The apparent 27r peri- 
odicity seems surprising, but occured also in the original 
calculation of a similar mode by Lieb [|| within a full 
quantum treatment of a ID Bose gas interacting via a 
(^-function potential. The Lieb mode was later rederived 
by a fairly accurate semiclassical approximation based 
on the elementary solitary wave of the ID classical GP 
model |,|. 

Lieb further argued that the corresponding Bogoliubov 
mode is no more elementary and thus proposed an in- 
triguing dual interpretation of the excitation spectrum. 
It should be noted that the dispersions of the two modes 
exhibit the same linear dependence at low momenta, 
E w c\Q\, where c is the speed of sound, but significant 
differences arise at finite momenta. The differences are 
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FIG. 6: Energy E in units of "/±(fiuJx) versus impulse Q in 
units of y±(h/aj_) for 7=1. The solid line corresponds to the 
complete sequence of solitary waves discussed in the text, and 
the dotted line to the 7=1 Bogoliubov dispersion of Fig. |^ 
adjusted to current units. Symbols I and II correspond to the 
two special cases of the solitary wave illustrated in Figs. ^ and 

I 

especially pronounced in the current calculation within 
a cylindrical trap. Specifically, let us assume an average 
linear density v — 0.02 atoms/A which leads to 7 = va= 1 
and 71 = va±_ = 10 2 . If we then adjust the 7=1 Bogoli- 
ubov dispersion of Fig. || to the units employed in Fig. 
the two dispersions are seen to diverge very quickly at 
the scale of Fig. |[ In other words, Bogoliubov and Lieb 
modes operate at rather different energy and momentum 
scales in a realistic trap. 

To summarize the preceding accurate calculation for 
7=1, the solitary wave is essentially quasi-lD in this 
weak-coupling region and its main features are indeed 
captured by an effective ID model EM. However, all 
cases of actual experimental interest f[7| are char- 

acterized by significantly larger values of the effective 
coupling where quasi-lD solitons are expected to be un- 
stable |l^, In particular, Ref. |l2| suggests that 
a critical coupling occurs in a cylindrical trap when 
n mEiX Uo / Hlu ± = 2.4, where n max is the maximum local 
particle density in the ground state of the trap. If we 
tentatively assume that the TF approximation (Q) can 
be trusted at the anticipated critical coupling, the above 
criticality condition reads 27 1 / 2 = 2.4 where 7 is the di- 
mensionless effective coupling constant defined in Eq. (jl]) . 
The root of this equation 7 = 7 C = 1-44 provides a critical 
coupling 7 C above which quasi-lD solitons are predicted 
to be unstable |3]J| . 

In fact, our numerical calculation suggests that the 



FIG. 7: Energy E in units of ^y±(hui±) versus impulse Q in 
units of j±(h/a±) for 7 = 10. The solid line corresponds to 
the fundamental sequence of vortex rings discussed in the text, 
and the dotted line to the auxiliary sequence that contains the 
black soliton (point IV). 



quasi- ID nature of the solitary wave is completely lost 
at a higher critical coupling, namely for 7 > 71 ps 3.9. 
The emerging new picture is clear at 7= 10 which is the 
special case described in our recent short communication 
. This case is reanalyzed and further extended in the 
continuation of the present paper. 

It is natural to begin again with the calculation of a 
static (v = 0) soliton obtained by using the input con- 
figuration ( |19[) with c\ = 0, C2 = 1 and practically any 
C3. We then increment the velocity to both positive and 
negative values in steps of <5v = ±0.01 to yield a sequence 
of solitary waves which now display two surprising fea- 
tures. First, a ringlike structure develops for 7= 10 that 
was not present at 7 = 1. Second, the above sequence 
exists only over the limited velocity range —V\ < v < Vi 
where v\ — 0.84 = 0.47 c and c = 1.77 is the speed of 
sound calculated in Sec. II for 7 = 10. The existence of a 
critical velocity v± also becomes apparent in the energy- 
momentum dispersion of the above sequence depicted by 
a dotted line in Fig. 0. This portion of the dispersion 
is symmetric around Q — n, where it achieves a maxi- 
mum, but remains open ended at two critical points that 
correspond to v = ±«i . 

It is thus not surprising that an independent sequence 
of solitary waves with lower energy exists for 7= 10. In- 
deed, we return to the input configuration ( J10| ) but now 
target a solution with velocity in the range v± < v < c. Af- 
ter some experimentation a solution is obtained for, say, 
v = 1.5 if we choose the trial parameters ci =0.2, C2 = 0.98 
and C3 = 3. Having thus obtained a specific solution for 
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FIG. 8: Solitary wave for 7= 10 and a representative value of 
the velocity v = c/2. We show the radial dependence of the 
real and the imaginary part of the wave function for various 
values of £ in steps of <5£ = 0.1. The corresponding results for 
negative £ are obtained through the parity relations (53). 



v = 1.5 the algorithm is iterated forward and backwards 
in steps of 5u = ±0.01 to obtain an entirely new sequence 
of solitary waves in the velocity range —v\ < v < c, and 
a corresponding sequence for — c < v < v% through the 
symmetry relation (pi]). Here v\ is the same critical ve- 
locity encountered in the preceding paragraph, as is also 
apparent in the calculated energy-momentum dispersions 
which are depicted by solid lines in Fig. |?j and join the 
previously calculated dotted line through cusps that cor- 
respond to V = ±V\. 

Hence we turn to a description of the detailed nature 
of this new sequence of solitary waves. For values of the 
velocity near the speed of sound c, the calculated soli- 
ton appears again as a mild soundlike disturbance of the 
ground state. The dominant features of the solitary wave 
are pronounced as the velocity is decreased to lower val- 
ues and become reasonably apparent for v = c/2 that 
corresponds to point I in the dispersion of Fig.[7| The 
wave function is completely illustrated through its real 
and imaginary parts in Fig. ||. An important new fea- 
ture emerges by comparison with the corresponding case 
at 7 = 1 illustrated in frame I of Fig. ^. Both the real and 
the imaginary part at the center of the soliton (£ = 0) now 
vanish for a specific radius R — 2.8, thus a vortex ring is 
beginning to emerge. A partial but more transparent il- 
lustration is given in Fig. || where we depict the radial 
dependence of the local particle density n = \^\ 2 for var- 
ious values of £. Again it is clear that the density near 
the center of the soliton (£ = 0) vanishes on a ring with 
a relatively large radius R = 2.8 . The features of the 




5 £=° 



FIG. 9: Radial dependence of the local particle density 
n — |vp| 2 for 7 = 10, using the same conventions for the £ 
dependence as in Fig. ^| The four special cases considered 
correspond to the four representative points I, IL III, and IV 
along the energy-momentum dispersion of Fig. bj 
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FIG. 10: Contour levels of the local particle density n=|\£ f | 2 
for 7 = 10, in a plane that contains the z axis and cuts across 
the cylindrical trap. The complete 3D picture may be envis- 
aged by simple revolution around the z axis. Regions with 
high particle density are bright while regions with zero den- 
sity are black. The four special cases considered are the same 
as in Fig. 0. 
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vortex ring become completely apparent, and its radius 
is tightened, as we proceed to smaller values of the ve- 
locity. A notable special case is the static (v = 0) vortex 
ring with radius R=1.8 illustrated in frame II of Fig. ||, 
which is far from being a black soliton. The correspond- 
ing point II in Fig. [7] is thus a new local maximum of 
the energy-momentum dispersion, which is clearly dis- 
tinguished from the local maximum at point IV that cor- 
responds to the static black soliton discussed earlier in 
the text. 

One would think that pushing the velocity v to neg- 
ative values would somehow retrace the calculated se- 
quence of vortex rings backwards. In fact, our algorithm 
continues to converge to vortex rings of smaller radii un- 
til the critical velocity —V\ is encountered where the ring 
achieves its minimum radius i? m in = 0.8 and ceases to ex- 
ist for smaller values of v. The terminal state at v — —V\ is 
illustrated in frame III of Fig. ||. We have thus described 
a sequence of solitary waves that consists of bonafide 3D 
vortex rings and does not contain a black soliton. The 
corresponding branch in the energy-momentum disper- 
sion of Fig. ^ is labeled by points I, II, and III that stand 
for the special cases i> = c/2,0, and —v%. As mentioned 
already, an equivalent sequence of solitary waves exists 
in the range — c < v < v\ and leads to a dispersion curve 
in Fig. |?j that is mirror symmetric to the branch (I, II, III) 
around Q = n. 

To complete the description for 7 = 10 we must briefly 
return to the auxiliary sequence of solitary waves associ- 
ated with the portion of the dispersion that is depicted 
by a dotted line in Fig. [?]. As one moves from point III to 
point IV, the ringlike structure is more or less preserved 
at constant radius R — R min = 0.8. Nevertheless, the 
detailed features of the vortex ring are tamed at small 
velocities and completely disappear for v — to yield a 
black soliton at point IV. 

We thus essentially conclude our description of soli- 
tary waves for 7 = 10 by schematically summarizing our 
main results in Fig. [l(j. Yet some of the elements of 
the preceding discussion are sufficiently surprising to de- 
serve closer attention. For example, simple inspection 
of Fig. fj] reveals that the group velocity becomes nega- 
tive in the region (II, III) or, equivalently, the impulse is 
opposite to the group velocity. This rotonlike behavior 
is consistent with the Onsager-Feynman view of a ro- 
ton as the ghost of a vanished vortex ring |52| because 
the calculated radius of the vortex ring is monotonically 
decreasing along the fundamental (I, II, III) sequence. A 
full-scale roton would develop if the terminal point III 
were an inflection point beyond which the group velocity 
begins to rise again. Actually, this is exactly what hap- 
pens as one moves away from point III along the upper 
branch in Fig. 0, but this "roton" portion of the disper- 
sion now appears in a strange location by comparison to 
the usual situation in liquid helium On the other 

hand, the black soliton at the stationary point IV is in- 
deed the ghost of a vanished vortex ring, as explained in 
the preceding paragraph. 




FIG. 11: Energy E in units of ~/±(hu)±) versus impulse Q 
in units of y±(h/a±) for 7 = 20. The solid lines correspond 
to the fundamental single-ring sequence, the dashed lines to 
the double-ring sequence, and the dotted line to the auxiliary 
sequence that contains a black soliton (point V). 



It is also interesting to question how the picture just 
described evolves with increasing values of the dimen- 
sionless coupling constant 7 which is the only parameter 
that enters the rationalized GP equation. Our numerical 
calculations have revealed yet another critical coupling 
72 ~ 12, in the sense that new flavor arises for 7 > 72- 
The structure of the solitary waves in this new regime 
becomes sufficiently clear for 7 = 20 and is best sum- 
marized by the calculated energy-momentum dispersion 
shown in Fig. 1 1 . Apart from mirror symmetry, the dis- 
persion now exhibits two cusps that correspond to two 
critical velocities v\ — 1.35 = 0.64c and V2 — 0.48 = 0.23 c, 
where c = 2.1 is the speed of sound calculated for 7 = 20 
as in Sec. II. 

The nature of the solitary waves associated with the 
various branches in the dispersion of Fig. [ll]is very briefly 
described with the aid of Fig. [l^. Thus we consider 
the sequence of five characteristic points (I,II,...,V) that 
roughly cover half of the dispersion, the other half be- 
ing obtained by the mirror symmetry (|2l|). The lowest 
branch (I, II) corresponds to single vortex rings with ve- 
locities in the range c > v > — v\, as in the case 7 = 10. 
Again the ring achieves its minimum radius at the critical 
velocity v — —v\ (point II). The new element for 7 = 20 is 
the intermediate branch (II, III, IV) that corresponds to 
double rings with velocities in the range v 2 > v > —v\. 
The second ring is first created at the flanks of the trap 
and comes closest to the original ring at the new critical 
velocity v=v% (point IV). This double-ring configuration 
is more or less preserved along the upper branch (IV, V) , 
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FIG. 12: Radial dependence of the local particle density 
n = l^l 2 for 7 = 20. The four special cases considered corre- 
spond to points I, II, III, and IV along the energy-momentum 
dispersion of Fig. till. Point V in the above dispersion in not 
illustrated here because it corresponds to a black soliton sim- 
ilar to that shown earlier in frame IV of Fig. [)| 



with velocities in the range vi > v > 0, but gradually fades 
away to become a black soliton at v = (point V). 

While the numerical calculation becomes increasingly 
more difficult for larger values of 7, it is clear that a 
sequence of critical couplings 71 , 72 , . . . exists and leads 
to a hierarchy of axisymmetric vortex rings. The single- 
ring solution associated with the lowest branch of the 
spectrum is a robust feature for all 7 > 71 ~ 3.9 and 
is likely stable to all perturbations. But it is possible 
that the multiple-ring configurations associated with the 
higher branches are unstable to non-axisymmctric per- 
turbations. In this respect, one should recall that the 
solitary wave that corresponds to the upper branch in 
the original calculation of Jones and Roberts jig] within 
the homogeneous GP model was later argued to be un- 
stable nf. 

However, we should emphasize that the vortex rings 
constructed here differ significantly from the Jones- 
Roberts (JR) vortex ring that provided the basic mo- 
tivation for the present work. As with ordinary smoke 
rings in fluid dynamics, the JR ring can never be static 
thanks to a virial relation of the type (|l6| ) that prevents 
finite-energy solutions with v = in the homogeneous GP 
model. As a result, the radius of the vortex ring grows 
to infinity at low velocity. This picture is completely re- 
arranged in a cylindrical trap because the occurrence of 
slow vortex rings with large radius is restricted by the 
boundaries of the trap. Instead, vortex rings are pre- 
dicted to nucleate at the flanks of the trap as soundlikc 



pulses with high velocity approaching the speed of sound 
c, and their radius actually decreases with decreasing ve- 
locity. In particular, it is now possible to obtain static 
(v = 0) vortex rings of finite radius that arc no longer 
contradicted by the virial relation ([h]) . The structure of 
the energy-momentum dispersions calculated throughout 
the present paper clearly reflects the substantial restruc- 
turing of vortex rings within a cylindrical trap. 

It is then natural to question whether or not there 
exists a limit in which the JR vortex ring is recovered. 
One should expect that this may happen when the bulk 
healing length L=(na) -1 / 2 is significantly smaller than 
the transverse oscillator length ax- This limit is trans- 
lated into large values of 7 = va which is the only di- 
mensionless parameter that enters the rationalized GP 
equation. Now, our discussion earlier in this section sug- 
gests an increasingly complicated hierarchical structure 
in the strong-coupling limit rather that a simple JR soli- 
ton. A logical conclusion is that a JR vortex ring some- 
how created within the bulk will sooner or later sense the 
boundaries of the cylindrical trap. It will thus either di- 
rectly dissipate into sound waves, or reorganize itself to 
conform with one or more of the presently calculated vor- 
tex rings possibly after ejecting some amount of radiation 
in the form of sound waves. 

The preceding remarks indicate a certain non- 
uniformity that is inherent in the approximation of the 
cigar-shaped trap by an infinitely-long cylindrical trap. 
The same phenomenon is also apparent in the calculation 
of linear modes in Sec. II. For instance, neither one of the 
two asymptotes for the speed of sound quoted in Eqs. ( p~l] ) 
and ([12|) approaches the well-known Bogoliubov speed 
in a homogeneous Bose gas jn]]. But a non-uniformity 
of this type is not a reason to doubt that a sufficiently 
elongated trap can be approximated by an infinitely-long 
cylindrical trap. 



IV. CONCLUSION 

We have thus presented a complete 3D calculation of 
solitary waves in a cylindrical Bose-Einstein condensate. 
In all cases considered there exists a nontrivial phase 
difference 6(j> that is reminiscent of strictly ID solitons 
|^|, D and is important for their experimental production 
through phase imprinting || ^ . 

Nevertheless, the detailed structure of the solitary 
waves depends crucially on the strength of the dimen- 
sionless effective coupling constant 7. Quasi-ID solitons 
occur only in the weak-coupling region 7 < 71 sa 3.9 where 
some of our accurate numerical results could be approx- 
imated through effective ID models But a suffi- 
ciently strong coupling or density is necessary in order to 
pronounce the special features of a condensate. It is thus 
not surprising that the effective coupling in experiments 
performed so far lies in the region 7 > 71 where the na- 
ture of the theoretically predicted solitary waves changes 
drastically. 
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For 7 > 71 solitary waves are still characterized by a 
nontrivial phase difference 5<f> between the two ends of 
the trap but are otherwise 3D vortex rings. This finding 
is consistent with a recent experiment |17j and a cor- 
responding theoretical analysis Jl3| in finite traps. The 
main mathematical advantage of the approximation of a 
sufficiently elongated trap by an infinitely-long cylindri- 
cal trap is that vortex rings can then be calculated in a 
steady state propagating with a constant velocity v. It is 
thus possible to carry out a detailed study of the soliton 
profile as a function of the effective coupling constant 7 
and the velocity v, as is done in the present paper. 

An interesting by-product of the above idealization 
is that a soliton is characterized by a definite energy- 
momentum dispersion. The calculated dispersion is 
found to be the direct analog of the Lieb mode || in the 
weak-coupling region and acquires interesting rotonlike 
features for stronger couplings. Perhaps such a disper- 
sion can be measured by a combination of phase imprint- 



ing ||, m with Bragg spectroscopy recently employed for 
the detection of the usual Bogoliubov mode |34[ pq] . In 
this respect, one should keep in mind that Bogoliubov 
and Lieb modes operate at different energy and momen- 
tum scales in a realistic trap. This fact becomes evident 
by the different sets of physical units employed for the 
Bogoliubov mode in Fig. ^| and the Lieb mode in, say, 
Fig. ||. The difference is accounted for by the second di- 
mensionless coupling 7j_ = va^ in Eq. (|]J) which is much 
stronger than j — i/a because 7j_/7 = a±/a~ 10 2 . 
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